Setup R notebook Environment

knitr::opts_chunk$set(echo = TRUE)
Sys.setlocale("LC_CTYPE", "en_US.UTF-8")
[1] ""

1 Package load for Project

packages = c("sp","sf","tidyverse","tmap","rjson", "rgdal", "geojsonio", "GISTools",'spatialEco', "raster","dplyr","spdep") 
for (p in packages){
  if(!require(p, character.only = T)){
    install.packages(p) 
  } 
  library(p,character.only = T) 
}
Loading required package: sf
Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3
Loading required package: tidyverse
-- Attaching packages --------------------------------------- tidyverse 1.2.1 --
v ggplot2 3.1.0     v purrr   0.2.5
v tibble  1.4.2     v dplyr   0.7.8
v tidyr   0.8.2     v stringr 1.3.1
v readr   1.3.1     v forcats 0.3.0
-- Conflicts ------------------------------------------ tidyverse_conflicts() --
x dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()
Loading required package: tmap
Loading required package: rjson
Loading required package: rgdal
rgdal: version: 1.3-6, (SVN revision 773)
 Geospatial Data Abstraction Library extensions to R successfully loaded
 Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
 Path to GDAL shared files: C:/Users/terence.tan.2015/Documents/R/win-library/3.5/rgdal/gdal
 GDAL binary built with GEOS: TRUE 
 Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
 Path to PROJ.4 shared files: C:/Users/terence.tan.2015/Documents/R/win-library/3.5/rgdal/proj
 Linking to sp version: 1.3-1 
Loading required package: geojsonio

Attaching package: <U+393C><U+3E31>geojsonio<U+393C><U+3E32>

The following object is masked from <U+393C><U+3E31>package:base<U+393C><U+3E32>:

    pretty

Loading required package: GISTools
Loading required package: maptools
Checking rgeos availability: TRUE
Loading required package: RColorBrewer
Loading required package: MASS

Attaching package: <U+393C><U+3E31>MASS<U+393C><U+3E32>

The following object is masked from <U+393C><U+3E31>package:dplyr<U+393C><U+3E32>:

    select

Loading required package: rgeos
rgeos version: 0.4-2, (SVN revision 581)
 GEOS runtime version: 3.6.1-CAPI-1.10.1 
 Linking to sp version: 1.3-1 
 Polygon checking: TRUE 

Loading required package: spatialEco
Loading required package: raster

Attaching package: <U+393C><U+3E31>raster<U+393C><U+3E32>

The following objects are masked from <U+393C><U+3E31>package:MASS<U+393C><U+3E32>:

    area, select

The following object is masked from <U+393C><U+3E31>package:dplyr<U+393C><U+3E32>:

    select

The following object is masked from <U+393C><U+3E31>package:tidyr<U+393C><U+3E32>:

    extract

Loading required package: spdep
Loading required package: Matrix

Attaching package: <U+393C><U+3E31>Matrix<U+393C><U+3E32>

The following object is masked from <U+393C><U+3E31>package:tidyr<U+393C><U+3E32>:

    expand

Loading required package: spData
To access larger datasets in this package, install the spDataLarge package with:
`install.packages('spDataLarge', repos='https://nowosad.github.io/drat/', type='source')`

Attaching package: <U+393C><U+3E31>spData<U+393C><U+3E32>

The following object is masked _by_ <U+393C><U+3E31>.GlobalEnv<U+393C><U+3E32>:

    x

The following object is masked from <U+393C><U+3E31>package:spatialEco<U+393C><U+3E32>:

    elev

2 Importing Taiwan GDAM Subzone

Shapefile Object with Spatial data Import data for past 12 years

Requires translation of chinese character independetly (https://stat.ethz.ch/pipermail/r-sig-geo/2009-March/005323.html)

Taiwan follows EPSG 3826 (http://spatialreference.org/ref/epsg/twd97-tm2-zone-121/)

taiwan_ts_map_sp <- readOGR(dsn = "data/TAIWAN_TOWNSHIP", layer = "TOWN_MOI_1071226", stringsAsFactors=TRUE)
OGR data source with driver: ESRI Shapefile 
Source: "C:\Users\terence.tan.2015\Desktop\dangy\data\TAIWAN_TOWNSHIP", layer: "TOWN_MOI_1071226"
with 368 features
It has 7 fields
df_dengue <- jsonlite::fromJSON("data/dengue_case.json") %>% 
  filter(as.Date(Onset_day) >= "2016-01-01" & as.Date(Onset_day)< "2017-01-01")
# Check for NA values for coordinates
 sum(is.na(df_dengue$Minimum_statistical_area_center_point_X))
[1] 0
 sum(is.na(df_dengue$Minimum_statistical_area_center_point_Y))
[1] 0
 df_dengue<- df_dengue[!is.na(df_dengue$Minimum_statistical_area_center_point_X),]
 df_dengue<- df_dengue[!(df_dengue$Minimum_statistical_area_center_point_X == 'None'),]
 
 # Type conversion
 df_dengue[, c(10,11,19,23,24)] <- sapply(df_dengue[, c(10,11,19,23,24)], as.numeric)
 
  # Transform into SF object
 sf_dengue <- st_as_sf(df_dengue, 
                       coords = c("Minimum_statistical_area_center_point_X",
                                  "Minimum_statistical_area_center_point_Y"),
                        crs = "+init=epsg:3826 +proj=longlat +ellps=WGS84 +no_defs",na.fail=FALSE)
sf_dengue <- na.omit(sf_dengue)
sf_dengue <- as(sf_dengue, 'Spatial')

Map name dict to map

#plot(taiwan_ts_map_sp, border="lightgrey")
#plot(transaction.boundary.poly, col="red",add=TRUE)
tmap_mode("plot")
tmap mode set to plotting
map <- tm_shape(taiwan_ts_map_sp)+
      tm_polygons('TOWNENG')+
       tm_borders(col = "grey40",alpha=0.5)+
  tm_view(alpha = 1)+
  tm_layout(basemaps = c('OpenStreetMap'))+
  tm_shape(sf_dengue)+
  tm_bubbles(col ="red",size = 0.1)
As of version 2.0, basemaps and basemaps.alpha have to be called from tm_view
#map 
#tmap_mode("view")
#map 

3 FEATURE 5 & 6

projection(taiwan_ts_map_sp) <- sf_dengue@proj4string
taiwan_ts_map_sp@data$poly.ids <- 1:nrow(taiwan_ts_map_sp)
#join poly and points
pts.poly <- point.in.poly(sf_dengue, taiwan_ts_map_sp)
although coordinates are longitude/latitude, st_intersects assumes that they are planar
#assign unique id per poly, full info
#descendng order of town with most
selected_aggregated_temp <-pts.poly@data %>% 
  dplyr::select(poly.ids, TOWNNAME, TOWNENG, Onset_day, Case_study_date, Notification_day) %>% 
  mutate(MONTH = months(as.Date(pts.poly@data$Onset_day))) %>%
  group_by(poly.ids,TOWNNAME, TOWNENG,MONTH) %>% 
  summarise(total_cases = n())
total_case_district <- selected_aggregated_temp %>%
   group_by(poly.ids) %>% 
    summarise(total_cases = n())
total_case_district <- merge(taiwan_ts_map_sp, total_case_district, by=c("poly.ids")) 
total_case_district@data <- total_case_district@data %>% 
  mutate(total_cases = ifelse(is.na(total_cases),0,total_cases))
tmap_mode("view")
tmap mode set to interactive viewing
tm_shape(total_case_district) +
  tm_fill(col="total_cases",id="TOWNENG",style="jenks")

4 Spatial Autocorrelation

taiwan_q <- poly2nb(taiwan_ts_map_sp, queen=TRUE)
plot(taiwan_ts_map_sp, border="lightgrey") 

plot(taiwan_q, coordinates(taiwan_q), pch = 19, cex = 0.6, add = TRUE, col= "red") 
Error in (function (classes, fdef, mtable)  : 
  unable to find an inherited method for function ‘coordinates’ for signature ‘"nb"’
LS0tDQp0aXRsZTogIlBST0pFQ1QiDQphdXRob3I6ICJKZXJyeSBUb2h2YW4sIEFuZyBLYWggRW5nLCBUZXJlbmNlIFRhbiwgSmVycnkgVG9odmFuIg0KZGF0ZTogIjEvMTQvMjAxOSINCm91dHB1dDogDQogIGh0bWxfbm90ZWJvb2s6DQogICAgbnVtYmVyX3NlY3Rpb246IHllcw0KICAgIHRoZW1lOiBmbGF0bHkNCiAgICB0b2M6IHllcw0KICAgIHRvY19mbG9hdDogeWVzDQogIHBkZl9kb2N1bWVudDoNCiAgICB0b2M6IHllcw0KICB3b3JkX2RvY3VtZW50Og0KICAgIHRvYzogeWVzDQotLS0NCg0KU2V0dXAgUiBub3RlYm9vayBFbnZpcm9ubWVudA0KYGBge3Igc2V0dXAsIGluY2x1ZGU9VFJVRSwgZXZhbD1UUlVFLG1lc3NhZ2U9RkFMU0UsIHdhcm5pbmc9RkFMU0V9DQprbml0cjo6b3B0c19jaHVuayRzZXQoZWNobyA9IFRSVUUpDQpTeXMuc2V0bG9jYWxlKCJMQ19DVFlQRSIsICJlbl9VUy5VVEYtOCIpDQpgYGANCg0KI1BhY2thZ2UgbG9hZCBmb3IgUHJvamVjdA0KYGBge3J9DQpwYWNrYWdlcyA9IGMoInNwIiwic2YiLCJ0aWR5dmVyc2UiLCJ0bWFwIiwicmpzb24iLCAicmdkYWwiLCAiZ2VvanNvbmlvIiwgIkdJU1Rvb2xzIiwnc3BhdGlhbEVjbycsICJyYXN0ZXIiLCJkcGx5ciIsInNwZGVwIikgDQpmb3IgKHAgaW4gcGFja2FnZXMpew0KICBpZighcmVxdWlyZShwLCBjaGFyYWN0ZXIub25seSA9IFQpKXsNCiAgICBpbnN0YWxsLnBhY2thZ2VzKHApIA0KICB9IA0KICBsaWJyYXJ5KHAsY2hhcmFjdGVyLm9ubHkgPSBUKSANCn0NCmBgYA0KDQojSW1wb3J0aW5nIFRhaXdhbiBHREFNIFN1YnpvbmUNClNoYXBlZmlsZSBPYmplY3Qgd2l0aCBTcGF0aWFsIGRhdGENCkltcG9ydCBkYXRhIGZvciBwYXN0IDEyIHllYXJzDQoNClJlcXVpcmVzIHRyYW5zbGF0aW9uIG9mIGNoaW5lc2UgY2hhcmFjdGVyIGluZGVwZW5kZXRseSAoaHR0cHM6Ly9zdGF0LmV0aHouY2gvcGlwZXJtYWlsL3Itc2lnLWdlby8yMDA5LU1hcmNoLzAwNTMyMy5odG1sKQ0KDQpUYWl3YW4gZm9sbG93cyBFUFNHIDM4MjYgKGh0dHA6Ly9zcGF0aWFscmVmZXJlbmNlLm9yZy9yZWYvZXBzZy90d2Q5Ny10bTItem9uZS0xMjEvKQ0KYGBge3J9DQp0YWl3YW5fdHNfbWFwX3NwIDwtIHJlYWRPR1IoZHNuID0gImRhdGEvVEFJV0FOX1RPV05TSElQIiwgbGF5ZXIgPSAiVE9XTl9NT0lfMTA3MTIyNiIsIHN0cmluZ3NBc0ZhY3RvcnM9VFJVRSkNCmRmX2Rlbmd1ZSA8LSBqc29ubGl0ZTo6ZnJvbUpTT04oImRhdGEvZGVuZ3VlX2Nhc2UuanNvbiIpICU+JSANCiAgZmlsdGVyKGFzLkRhdGUoT25zZXRfZGF5KSA+PSAiMjAxNi0wMS0wMSIgJiBhcy5EYXRlKE9uc2V0X2RheSk8ICIyMDE3LTAxLTAxIikNCg0KIyBDaGVjayBmb3IgTkEgdmFsdWVzIGZvciBjb29yZGluYXRlcw0KIHN1bShpcy5uYShkZl9kZW5ndWUkTWluaW11bV9zdGF0aXN0aWNhbF9hcmVhX2NlbnRlcl9wb2ludF9YKSkNCiBzdW0oaXMubmEoZGZfZGVuZ3VlJE1pbmltdW1fc3RhdGlzdGljYWxfYXJlYV9jZW50ZXJfcG9pbnRfWSkpDQogZGZfZGVuZ3VlPC0gZGZfZGVuZ3VlWyFpcy5uYShkZl9kZW5ndWUkTWluaW11bV9zdGF0aXN0aWNhbF9hcmVhX2NlbnRlcl9wb2ludF9YKSxdDQogZGZfZGVuZ3VlPC0gZGZfZGVuZ3VlWyEoZGZfZGVuZ3VlJE1pbmltdW1fc3RhdGlzdGljYWxfYXJlYV9jZW50ZXJfcG9pbnRfWCA9PSAnTm9uZScpLF0NCiANCiAjIFR5cGUgY29udmVyc2lvbg0KIGRmX2Rlbmd1ZVssIGMoMTAsMTEsMTksMjMsMjQpXSA8LSBzYXBwbHkoZGZfZGVuZ3VlWywgYygxMCwxMSwxOSwyMywyNCldLCBhcy5udW1lcmljKQ0KIA0KICAjIFRyYW5zZm9ybSBpbnRvIFNGIG9iamVjdA0KIHNmX2Rlbmd1ZSA8LSBzdF9hc19zZihkZl9kZW5ndWUsIA0KICAgICAgICAgICAgICAgICAgICAgICBjb29yZHMgPSBjKCJNaW5pbXVtX3N0YXRpc3RpY2FsX2FyZWFfY2VudGVyX3BvaW50X1giLA0KICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICJNaW5pbXVtX3N0YXRpc3RpY2FsX2FyZWFfY2VudGVyX3BvaW50X1kiKSwNCiAgICAgICAgICAgICAgICAgICAgICAgIGNycyA9ICIraW5pdD1lcHNnOjM4MjYgK3Byb2o9bG9uZ2xhdCArZWxscHM9V0dTODQgK25vX2RlZnMiLG5hLmZhaWw9RkFMU0UpDQpzZl9kZW5ndWUgPC0gbmEub21pdChzZl9kZW5ndWUpDQpzZl9kZW5ndWUgPC0gYXMoc2ZfZGVuZ3VlLCAnU3BhdGlhbCcpDQpgYGANCg0KTWFwIG5hbWUgZGljdCB0byBtYXANCmBgYHtyfQ0KDQojcGxvdCh0YWl3YW5fdHNfbWFwX3NwLCBib3JkZXI9ImxpZ2h0Z3JleSIpDQojcGxvdCh0cmFuc2FjdGlvbi5ib3VuZGFyeS5wb2x5LCBjb2w9InJlZCIsYWRkPVRSVUUpDQoNCnRtYXBfbW9kZSgicGxvdCIpDQptYXAgPC0gdG1fc2hhcGUodGFpd2FuX3RzX21hcF9zcCkrDQogICAgICB0bV9wb2x5Z29ucygnVE9XTkVORycpKw0KICAgICAgIHRtX2JvcmRlcnMoY29sID0gImdyZXk0MCIsYWxwaGE9MC41KSsNCiAgdG1fdmlldyhhbHBoYSA9IDEpKw0KICB0bV9sYXlvdXQoYmFzZW1hcHMgPSBjKCdPcGVuU3RyZWV0TWFwJykpKw0KICB0bV9zaGFwZShzZl9kZW5ndWUpKw0KICB0bV9idWJibGVzKGNvbCA9InJlZCIsc2l6ZSA9IDAuMSkNCg0KI21hcCANCmBgYA0KDQoNCmBgYHtyfQ0KI3RtYXBfbW9kZSgidmlldyIpDQojbWFwIA0KDQpgYGANCg0KI0ZFQVRVUkUgNSAmIDYNCmBgYHtyfQ0KcHJvamVjdGlvbih0YWl3YW5fdHNfbWFwX3NwKSA8LSBzZl9kZW5ndWVAcHJvajRzdHJpbmcNCnRhaXdhbl90c19tYXBfc3BAZGF0YSRwb2x5LmlkcyA8LSAxOm5yb3codGFpd2FuX3RzX21hcF9zcCkNCg0KI2pvaW4gcG9seSBhbmQgcG9pbnRzDQpwdHMucG9seSA8LSBwb2ludC5pbi5wb2x5KHNmX2Rlbmd1ZSwgdGFpd2FuX3RzX21hcF9zcCkNCiNhc3NpZ24gdW5pcXVlIGlkIHBlciBwb2x5LCBmdWxsIGluZm8NCmBgYA0KDQpgYGB7cn0NCiNkZXNjZW5kbmcgb3JkZXIgb2YgdG93biB3aXRoIG1vc3QNCnNlbGVjdGVkX2FnZ3JlZ2F0ZWRfdGVtcCA8LXB0cy5wb2x5QGRhdGEgJT4lIA0KICBkcGx5cjo6c2VsZWN0KHBvbHkuaWRzLCBUT1dOTkFNRSwgVE9XTkVORywgT25zZXRfZGF5LCBDYXNlX3N0dWR5X2RhdGUsIE5vdGlmaWNhdGlvbl9kYXkpICU+JSANCiAgbXV0YXRlKE1PTlRIID0gbW9udGhzKGFzLkRhdGUocHRzLnBvbHlAZGF0YSRPbnNldF9kYXkpKSkgJT4lDQogIGdyb3VwX2J5KHBvbHkuaWRzLFRPV05OQU1FLCBUT1dORU5HLE1PTlRIKSAlPiUgDQogIHN1bW1hcmlzZSh0b3RhbF9jYXNlcyA9IG4oKSkNCg0KdG90YWxfY2FzZV9kaXN0cmljdCA8LSBzZWxlY3RlZF9hZ2dyZWdhdGVkX3RlbXAgJT4lDQogICBncm91cF9ieShwb2x5LmlkcykgJT4lIA0KICAgIHN1bW1hcmlzZSh0b3RhbF9jYXNlcyA9IG4oKSkNCg0KdG90YWxfY2FzZV9kaXN0cmljdCA8LSBtZXJnZSh0YWl3YW5fdHNfbWFwX3NwLCB0b3RhbF9jYXNlX2Rpc3RyaWN0LCBieT1jKCJwb2x5LmlkcyIpKSANCnRvdGFsX2Nhc2VfZGlzdHJpY3RAZGF0YSA8LSB0b3RhbF9jYXNlX2Rpc3RyaWN0QGRhdGEgJT4lIA0KICBtdXRhdGUodG90YWxfY2FzZXMgPSBpZmVsc2UoaXMubmEodG90YWxfY2FzZXMpLDAsdG90YWxfY2FzZXMpKQ0KDQp0bWFwX21vZGUoInZpZXciKQ0KdG1fc2hhcGUodG90YWxfY2FzZV9kaXN0cmljdCkgKw0KICB0bV9maWxsKGNvbD0idG90YWxfY2FzZXMiLGlkPSJUT1dORU5HIixzdHlsZT0iamVua3MiKQ0KYGBgDQoNCg0KDQoNCiNTcGF0aWFsIEF1dG9jb3JyZWxhdGlvbg0KYGBge3J9DQp0YWl3YW5fcSA8LSBwb2x5Mm5iKHRhaXdhbl90c19tYXBfc3AsIHF1ZWVuPVRSVUUpDQpwbG90KHRhaXdhbl90c19tYXBfc3AsIGJvcmRlcj0ibGlnaHRncmV5IikgDQpwbG90KHRhaXdhbl9xLCBjb29yZGluYXRlcyh0YWl3YW5fcSksIHBjaCA9IDE5LCBjZXggPSAwLjYsIGFkZCA9IFRSVUUsIGNvbD0gInJlZCIpIA0KYGBgDQoNCg==